Climate change induced habitat expansion of nutria (Myocastor coypus) in South Korea

The nutria, (Myocastor coypus), is a semiaquatic rodent native to the subtropical and temperate regions of South America. The species was introduced to South Korea for meat and fur production purposes and a wild population has become established. The species subsequently invaded aquatic ecosystems and destroyed aquatic vegetation and cultivated crops. Thus, it is essential to understand their current distribution and future range expansion for effective control and eradication strategies to reduce the risk of colonization into new regions. In this study, we used niche modeling procedure to identify potentially suitable habitats for M. coypus under current and future predicted climate change using the maximum entropy algorithm. We found that the main habitat area of M. coypus is expected to expand under a warming climate from ~ 4069 km2 in the southern and southeastern regions of South Korea, to the northern border of the country, with estimated ranges of 21,744 km2, 55,859 km2, and 64,937 km2 by 2030, 2050, and 2070, respectively. The findings of the present study assist in identifying the future distribution and potential dispersion routes of M. coypus in South Korea, which is important for informing the government regarding essential management actions plans at regional and local scales.


Results
Species occurrence survey. A species occurrence survey was conducted to support the "M. coypus Model performance and variable importance. Pearson's correlation test was performed, and based on the weak correlations (r < 0.70), five bioclimatic variables (e.g., annual mean temperature [Bio 1], maximum temperature of the warmest month [Bio 5], mean temperature of the wettest quarter [Bio 8], precipitation seasonality [ Bio 15], and precipitation of the wettest quarter [Bio 16]) and three environmental variables (i.e., altitude, distance from water, and land-cover change; Table S1) were selected and used for MaxEnt modeling. A jackknife test was performed to confirm the contribution of the included variables in the model (Fig. S1) and the results showed that the temperature-related variables Bio 5, Bio 8, and Bio 1 were the strongest bioclimatic predictors of the M. coypus distribution, explaining 29.23%, 21.05%, and 15.01% of the variance, respectively (Table 1). Similarly, both altitude and distance from water played relatively important roles, explaining 14. coypus habitats was modeled to show the distribution and the potential habitat area across all South Korean ADs under current and future climate change scenarios RCP 4.5, and RCP 8.5 22 in 2030, 2050, and 2070 ( Table 2). The current potential M. coypus habitat was identified in 12 provinces ( Fig. 2A), with a total area of ~ 4,069 km 2 , equivalent to 4.26% of the total country land area. The seven provinces concentrated along the southern and www.nature.com/scientificreports/  www.nature.com/scientificreports/ southeastern regions of the country, i.e., Gyeongsangnam (2,211 km 2 ), Jeollanam (506 km 2 ), Gyeongsangbuk (498 km 2 ), Busan (188 km 2 ), Daegu (213 km 2 ), Chungcheongbuk (160 km 2 ), and Jeollabuk (142 km 2 ), were predicted to be highly suitable for M. coypus habitation (Table 2). Suitability maps are independent of thresholds and commonly represent the probability of focal species to occur in a given geographical area. Therefore, we predicted the probability distribution maps of M. coypus and presented the current suitable habitats in Fig. 3A.  Seoul  0  385  563  577  399  595  600   Incheon  0  12  458  501  1  560  614   Gangwon  0  0  540  1,600 51  2,176 3,298   Gyeonggi  0  495  5,670 7,277 1,271 9,023 9,454   Sejong  5  4  160  334  26  397  437   Chungcheongbuk  160  32  673  1,951 93  3,329 4,430   Daejeon  0  39  321  469  61  479  526   Chungcheongnam  29  448  2,410 6,147 4A). Similarly, the mean suitability between the RCP 4.5 and RCP 8.5 for the years 2030, 2050, and 2070 was highest under RCP 8.5, suggesting that a warming climate will favor M. coypus habitation (Fig. 4B). Although the rates and trends of habitat expansion were inconsistent, suitable habitats were expected to increase across all provinces (except the cities of Ulsan and Busan) by 2030 under RCP 4.5. The areas with the highest suitability were predicted to be in Jeollanam, Gyeongsangnam, Gyeongsangbuk, Chungcheongnam, and Gyeonggi by 2030, 2050, and 2070 under both examined climate change scenarios (Table 2). However, the mean suitability value was estimated to be the highest in Gwangju City where it was predicted to range from P = 0.70-1.00 between 2030 and 2070 ( Table 3). The ADs located on the west coast (i.e., Jeollabuk, Chungcheongnam, and Gyeonggi) maintained higher rates of habitat expansion compared to those in the central region and along the east coast. Gangwon Province (northeastern region) and Jeju (southern region) had comparatively small areas of suitable habitat. In addition, high mountain regions (e.g., the Baekdu-daegan mountain range and Hall Mountain of Jeju [elevation, 1950 m]) also exhibited low habitat suitability (Fig. 3).

Vulnerability estimations of different administrative divisions.
We performed vulnerability estimation of different ADs based on the mean habitat suitability values of M. coypus in different provinces, which revealed that all the ADs were considered to have low vulnerability. However, in the future, vulnerability will be extended from the southern region to central and northern regions and high vulnerability with be exhibited in two ADs (e.g., Daejeon and Gwangju) by 2030, eight ADs (e.g., Seoul, Incheon, and Daejeon) by 2050, and 11 ADs (e.g., Seoul, Incheon, and Gyeonggi) by 2070 (Table 4).

Discussion
A 5-year field survey was performed to investigate the M. coypus distribution in South Korea, including identifying their presence in regions of the Nakdong River basin, Han River, and Jeju Island. The majority of habitat tracks were observed in the midstream and downstream regions of the Nakdong River basin, thus identifying this basin as the area that was the most populated by M. coypus. The midstream and downstream reaches of the Nakdong River are characterized by relatively slow water speeds upstream of several tributaries and wetlands 10   www.nature.com/scientificreports/ or Chungcheongbuk Province after 2014, implying significant population declines following the implementation of the eradication project. Globally, temperature is the most importnat factor for biological invasions 14 . In particular, M. coypus are warm-adapted species that are vulnerable to colder weather; therefore, they tend to propagate rapidly in areas with warm winters 23  As previously stated, the Nakdong River is currently the most suitable habitat for M. coypus in South Korea. Increasing temperatures facilitate invasive species growth, range expansion, and naturalization 27 , while simultaneously reducing the resilience of ecosystems to exotic species 14 . The model used here predicted that M. coypus will retain their current distribution and that additional habitat expansion will occur in the future, which is supported by the findings of Hilts et al. 14 and Pereira et al. 28 . Climatic factors generally have a dominant influence on species distribution at broad spatial scales 29 . The findings in the present study show that temperature-related variables are important determinants of climate-based M. coypus habitats (similar to the findings of Schertler et al. 13 . The results also predict habitat expansion from the southeastern region (particularly Gyeongsangnam, Jeollanam, Gyeongsangbuk, Busan City, and Daegu City) toward the central and northern regions of the country, which, in addition to Jeollabuk, Chungcheongnam, Chungcheongbuk, Gyeonggi, and Seoul City, were predicted to become highly vulnerable in future decades. These results support previous findings that show that increasing temperatures will expand the invasion risk northwards 30 , in addition to earlier projections indicating habitat expansion of invasive species with increasing temperatures 13,14,31 . Future expansion of the suitable M. coypus habitat is likely to occur gradually from south to north through the water channel connections of the five major rivers in South Korea. Northward expansion is expected to occur mainly from the west coast. The terrain of South Korea is characteristically low along the west coast and high along the east coast 31 . Consequently, most major rivers flow westward into the Yellow Sea and are relatively long, maintaining smooth slopes and wide basins. In contrast, the rivers flowing into the East Sea are relatively short and steeply sloped. The geographic terrain and river systems along the west coast, combined with increasing future temperatures, could favor the expansion of M. coypus habitat, supporting the future establishment of M. coypus populations. The derived model used in this study did not predict M. coypus habitats in high mountains or on sharp elevational gradients, such as the Baekdu-daegan mountain range (701 km), Mt. Halla on Jeju Island (1,950 m), the northern part of Gyeongsangbuk Province, and most of Gangwon Province (as similarly shown by Sheffels 23 and Pereira et al. 28 ). These areas could thus serve as a geographical barrier, blocking any future dispersion. In addition, high-elevational areas experience extremely cold weather in winter, maintaining snowpack for long time periods, thus creating inhospitable conditions for M. coypus 15 .
In addition to bioclimatic variables related to local-scale features, habitat covariates, such as a high proportion of freshwater and forested shrub wetlands (i.e., marshlands or swamps) close to other wetlands, appeared likely to support increasing M. coypus populations 14 . Freshwater swamps are relatively productive systems that enhance high floral diversity, including densely layered understories of herbaceous plants, shrubs, young trees, and overstory trees 32 , possibly generating foraging habitat and thermal protection during extremely cold winters. M. coypus tend to prefer various aquatic plants (e.g., Paspalum disthichum, Panicum tricholaenoides, Lemna minor, Spirodela polyrrhiza, Carex spp., and Schoenoplectus spp.), depending on their habitat and the season 33,34 . In South Korea, Hong et al. 19 reported M. coypus diets consisting of similar aquatic plants in the families of Poaceae, Cyperaceae, and Salviniacease, with some terrestrial plants belonging to the families of Chenopodiaceae, Asteraceae, and Amaranthaceae. Freshwater swamps, particularly those connected to groundwater, freshwater rivers, and lakes, could be appropriate M. coypus habitats considering the preference of this species for foraging near the edges of slow-speed water bodies, potentially maintaining connectivity between rivers and wetlands, facilitating  39 , potentially resulting in shortages of food resources for other wild herbivores inhabiting areas adjacent to M. coypus habitat (e.g., roe deer) 40 . M. coypus can weaken dams and irrigation systems by burrowing on riverbanks, thereby creating direct economic losses, as well as indirect losses connected to floods. Accordingly, the Ministry of Environment in South Korea conducted the M. coypus Eradication Project from 2014-2018, which was considered to have successfully reduced the range of M. coypus habitats in Jeju and Chungcheongbuk Provinces; however, major habitats persisted in the Nakdong River.
Based on the results of the present study, future suitable M. coypus habitats are likely to expand from the southeastern region to the northern international border. Currently, all the ADs are under the low vulnerable; however, by 2070, 11 ADs including Gwangju, Jeollabuk, Busan, Daegu, and Seoul, will be considered to have high vulnerability. The continuous M. coypus habitat expansion has already had adverse impacts on aquatic ecosystems and agriculture within South Korea, and the SDM-based prediction of M. coypus distribution provides a valuable tool for targeting areas with the highest risk of invasion. The results can assist with understanding M. coypus distribution and dispersion in South Korea, in addition to informing government and conservation agencies regarding management action plans and policies aimed at controlling or eradicating M. coypus at regional and local scales.

Methods
Study area. The study area included the mainland and all the islands of South Korea, with a total land mass of 100,033 km 2 (Fig. 5). The country is divided into 17 ADs (provinces and metro cities, referred to here as cities), which were utilized to assist in delineating M. coypus distribution localities and informing provincial governments when adopting future management actions. Geographically, the terrain of South Korea is mostly mountainous, dominating the northern and eastern parts of the country, whereas lowlands and flat plains occupy the southern and western regions. The Baekdu-daegan mountain range spans ~ 701 km, starting from Hyanro-bong in Goseong County of Gangwon Province and ending at Cheonwang-bong on Mt. Jiri in Jeollanam Province. The climate of South Korea is readily divisible into warm-temperate (southern coast and islands), temperate (northern and central), and cold-temperate (high mountains) 41 . Mean winter temperatures range between − 6 and 3 °C, whereas mean summer temperatures range from 23-26 °C 16 . Annual precipitation is between 1200-1500 mm, with a relatively higher rate of precipitation in the southern regions (including Jeju Island) compared www.nature.com/scientificreports/ to the central and northern regions 16 . The forests can be roughly categorized as deciduous broadleaf, temperate broadleaf, coniferous, subalpine, and alpine, with overall total species diversity of 41,483, including 5308 plant species, 1899 vertebrate species, and 22,612 invertebrate species (including 15,651 insect species) 42 .

M. coypus survey and species occurrence points.
A pilot survey was conducted to investigate the spatial distribution of M. coypus in South Korea based on earlier reports from the South Korean Ministry of Environment 43 and previous research 10 . The survey revealed the presence of M. coypus across all areas of the Nakdong River basin and its drainage system, in the surrounding lakes, in the downstream area of the Namhan River, and in the Jeju Island wetlands. Detailed surveys were then performed in the Nakdong River basin and its tributaries (127° 29ʹ-129° 19ʹ E and 35° 03ʹ-37° 13ʹ N), the Namhan River and Chungjuho artificial lake located downstream (127° 50ʹ-128° 20ʹ E and 36° 50ʹ-37° 05ʹ N), and across all the Jeju Island wetlands (126° 08ʹ-126° 58ʹ E and 33° 06ʹ-34° 00ʹN) (Fig. 5). Field surveys were performed from March 2014 to December 2018. During the surveys, all surveyors traveled on foot across river basins with an average walking speed of < 2 km·h −1 to monitor for footprints, living dens, feces, or make direct observations of individuals. Camera traps and M. coypus catching traps were installed to further confirm habitat and signs of presence. Geographic coordinates were recorded for all observations using a handheld geographical positioning system (GPS V; Garmin Inc., KS, USA). Habitat tracks were determined using the coordinate information within an area of 0.01 km 2 in ArcGIS 10.3 (Esri, CA, USA) 7 . GPS coordinate information for some representative photographs from the M. coypus survey can be seen in Fig. S2. The species occurrence points for each year are shown in Fig. 1. 44 , in addition to distance from water, altitude, and land-cover change (Table S2), were considered important for predicting the M. coypus distribution in South Korea. Climatic data, including precipitation and the monthly minimum or maximum temperatures, were obtained from the Korea Meteorological Administration (KMA; https:// www. kma. go. kr/ eng) to estimate the current climate and future climate change scenarios of South Korea. Two greenhouse gas emission scenarios were selected for analysis, i.e., RCP 4.5 and RCP 8.5, for 2030, 2050, and 2070. The moderate and maximum emission scenarios (RCP 4.5 and RCP 8.5, respectively) correlate to predicted global mean surface temperature increases of 1.4-1.8 °C and 2.0-3.7 °C by 2100, respectively. Among the numerous available global circulation models, the HadGEM3-RA is a regional atmospheric model produced by the Met Office Hadley Centre (https:// www. metoffi ce. gov. uk/), based on the atmospheric component of the current Earth System Model 45 . The KMA used HadGEM3-RA 45 to create a national climate change scenario for South Korea, revealing that it tended to model small-scale features more accurately than other global circulation models, such as HADGEM2-AO, which contains more complicated topographies, lengthy and uneven coastlines, and thousands of islands along the Korean Peninsula. Accordingly, the HadGEM3-RA global circulation model was used here to determine the impacts of the RCP 4.5 and RCP 8.5 scenarios, according to the R Package Dismo v. 1.3 46 . Similarly to Jeon et al. 47 , the current climate was determined by averaging climatic data recorded between 1950 and 2000, and future climates for 2030, 2050, and 2070 were estimated from predictions for 2026-2035, 2046-2055, and 2066-2075, respectively. We used a threshold of 0.7 for Pearson's correlation coefficient 48 for selecting bioclimatic and environmental (as described in Shin et al. 49 and Adhikari et al. 50 ). Seven landcover categories (i.e., urban, cropland, forest, grassland, wetland, barren, and water) were used in the model developed by Song et al. 51 , based on shared socioeconomic pathways using a scenario generator in which the future population and predicted urban area were adopted for a transition matrix containing the land-cover change trends of each class. The current and future land-cover change scenarios (shared socioeconomic pathways 1) developed by Song et al. 51 were obtained from the Korea Environment Institute MOTIVE (www. motive. kei. re. kr). Changes in the altitudinal range were quantified in the M. coypus distribution according to a digital elevation model. The variable distance from water included the distance from all kinds of rivers, lakes, and ponds, which were determined using the Euclidean distance option of the Spatial Analyst tool in ArcGIS Desktop 10.8. All the bioclimatic and environmental variables used in this study have a spatial resolution 0.01° (36 s), or ~ 1 km 2 , as explained by Fick and Hijmans 52 .

Environmental variables. Nineteen bioclimatic variables
Model development. SDM is a method used for predicting potential distribution of species throughout global space and time by using a correlation between the species' geographic occurrence and its surrounding environment 53 . SDMs typically estimate the ecological niche by statistically relating environmental variables directly to species occurrence. Among the various SDM algorithms, the MaxEnt is a widely used machine learning technique for habitat suitability, which has high predictive performance using only presence data 18 . In this study, we performed MaxEnt modeling using "Biomod2" Package v.3.5.1, selecting single model MAXENT. Phillips.2 54 to predict the current and future M. coypus distribution in South Korea. MaxEnt exhibits a high predictive performance using presence-only occurrence data 18,55 and is often the best option for invasive species because presence-absence survey data are usually not available, as was the case in the present study. As the Max-Ent model required background point data (e.g., pseudo-absence), we determined the background points of the study area using ArcGIS Desktop 10.8 (as suggested by Barbet-Massin et al. 56 ). Multiple species presence points within the same ~ 1 km 2 grid were eliminated, leaving only a single unique point per grid by adopting a spatially rarefy occurrence data tool in the ArcGIS SDM toolbox v.2.4 57 , to prevent overfitting and incorrect inflation of model performance owing to spatial autocorrelation 58 . The total species occurrence points for M. coypus were reduced to 2535 from 6180, and those points were used in the MaxEnt modeling. The five bioclimatic variables, distance from water, altitude, land-cover change, and the species occurrence points were used as the inputs for MaxEnt.
During modeling, the species occurrence data were randomly split into two parts, i.e., a training set and a test set, at a ratio of 3:1(as described in Adhikari et al. 59  Model evaluation and validation. The goodness-of-fit of the model was assessed, and a validation was performed with the AUC values of the receiver operating characteristic curves 61 and the TSS 62 . The training and test data sets were used to compute the AUC values. The AUC value measures the model capacity in a model that discriminates the observed occurrences from the background data using the training and test datasets 63 . AUC is independent of dataset size (prevalence) but is limited because it weights both commission and omission errors equally, avoids the actual probability values, and is dependent on geographical extent 64 . In particular, when expanding the geographical extent beyond the present range, the AUC values increase; thus, TSS was used as an alternative criterion for model performance validation. TSS corresponds to the modeled sum of sensitivity (omission error), which is the proportion of observed presences correctly predicted, and specificity (commission error), which is the proportion of observed absences correctly predicted 61 . AUC values vary from 0-1, whereby higher values suggest model superiority. Here, model performance was classified as poor (0.6-0.7), fair (0.7-0.8), good (0.8-0.9), or excellent (0.91-1), as suggested by Swets 65 . TSS is threshold-dependent and determines model performance by examining classification accuracy after selecting a threshold value. TSS values range between − 1 and + 1, where low values indicate an agreement that is no better than random chance, and high values represent the perfect alignment of observations and predictions 65 . In addition, a jackknife test was performed to evaluate the importance of each variable with respect to model performance 18,66 . The database and a brief conceptual flowchart of the methodology are presented in Fig. 6.

Prediction of potential habitat and new habitat expansion and estimation of vulnerabil-
ity. The MaxEnt model yielded both binary distribution maps and suitability maps of M. coypus under current and future climate change scenarios RCP 4.5 and RCP 8.5 for 2030, 2050, and 2070. The binary distribution maps were used to determine the suitable habitat areas across the 17 ADs of South Korea, and suitability maps were used to determine the mean and standard deviation of suitability values for each province. To calculate the distribution area, mean suitability values, and standard deviation of suitability values for each AD, the shape file representing 17 provinces of South Korea was overlain on the species distribution maps separately, and an extraction of the multiple values to points was performed using the zonal statistics of the spatial analyst tool in ArcGIS Desktop 10.8 31,67 . The vulnerability of each AD was determined and classified into three categories (low, moderate, or high) based on the mean suitability values from 0-0.33, 0.34-0.66, and 0.67-1.00, respectively. We used a linear scale to classify the vulnerability and the method of Adhikari et al. 59 with minor modifications. M. coypus habitat expansion was determined by differentiating between the current and future habitats using R package 'Raster' v3.5 68 as described by Jeon et al. 47 .
Ethics approval. Permission to collect species presence data are obtained from the Ministry of Environment, Republic of Korea.

Data availability
All data generated or analyzed during this study are included in this article.